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Abstract. A mechanism for fast magnetic reconnection in collisionless plasma is 
studied for understanding sawtooth collapse in tokamak discharges. Explosive growth 
of the tearing mode driven by electron inertia is analytically estimated by using an 
energy principle with a nonlinear displacement map. Decrease of the potential energy in 
the nonlinear regime (where the island width exceeds the electron skin depth) is found 
to be steeper than in the linear regime, resulting in accelerated reconnection. Release 
of free energy by such ideal fluid motion leads to unsteady and strong convective flow, 
which is not deterred by the small dissipation effects in high-temperature tokamak 
plasmas. Direct numerical simulation in slab geometry substantiates the theoretical 
prediction of the nonlinear growth. 
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1. Introduction 

Sawtooth collapse in tokamak plasmas has been a puzzling phenomenon for decades. 
Although the m = 1 kink-tearing mode is essential for the onset of this dynamics, 
Kadomtsev's full reconnection model [1] and the nonlinear growth of the resistive m = 1 
mode [2] (both based on resistive magnetohydrodynamic (MHD) theory) fail to explain 
the short collapse times (~ 100/is) as well as the partial reconnections observed in 
experiments [3l SI E] . Since resistivity is small in high-temperature tokamaks, two-fluid 
effects are expected to play an important role for triggering fast (or explosive) magnetic 
reconnection as in solar flares and magnetospheric substorms. 

In earlier works [6l [7], the linear growth rate of the kink-tearing mode in the 
collisionless regime has been analyzed extensively by using asymptotic matching, which 
shows an enhancement of the growth rate due to two-fluid effects, even in the absence 
of resistivity. Furthermore, direct numerical simulations [H [9l |10] of two-fluid models 
show acceleration of reconnection in the nonlinear phase, even though realistic two-fluid 
simulation of high-temperature tokamaks is still a computationally demanding task 
(especially when the resistive layer width is smaller than the electron skin depth de ~ 
1mm). These simulation studies, as a rule, indicate explosive tendencies of collisionless 
reconnection. 

However, theoretical understanding of such explosive phenomena is not yet 
established due to the lack of analytical development. In the neighborhood of the 
boundary (or reconnecting) layer, a perturbative approach breaks down at an early 
nonlinear phase and, consequently, asymptotic matching requires a fully nonlinear 
inner solution [11]. Moreover, in contrast to the quasi-equilibrium analysis developed 
for resistive reconnection [21 [12], the explosive process of collisionless reconnection 
should be a nonequilibrium problem, in which inertia is not negligible in the force 
balance and hence leads to acceleration of flow. Thus, the convenient assumption of 
steady reconnection is no longer appropriate. Recent theories [T3| [HI [15] emphasize 
the Hamiltonian nature of two-fluid models and try to gain deeper understanding of 
collisionless reconnection in the ideal limit. 

The purpose of the present work is to predict the explosive growth of the kink- 
tearing mode analytically by developing a new nonlinear variational technique that 
is based on a generalization of the MHD energy principle [T6| [T7] (for generalization 
see [IE])- For simplicity, we concentrate on the effect of electron inertia, which is 
an attractive mechanism for triggering fast reconnection in tokamaks; estimates of 
the reconnection rate are favorable [19], nonlinear acceleration is possible [9], and 
even the more mysterious partial reconnection may be explained by an inertia-driven 
collapse model [201 [21]. While we address the same problem as that of Ref. [9], our 
estimated nonlinear growth is quantitatively different from that of this reference, and 
our result is confirmed by direct numerical simulation. This advance in nonlinear theory 
is indispensable for clarifying the acceleration mechanism of collisionless reconnection. 

The present paper is organized as follows. In Sec. [21 we invoke a conventional 2D 
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slab model for electron inertia-driven reconnection and then construct its Lagrangian 
in terms of the fluid flow map (as in the ideal MHD theory [22| [23l [21] )• In Sec. [3l 
we obtain the linear growth rate of the inertial tearing mode in the large-A' regime 
(corresponding to the m = 1 kink-tearing mode in tokamaks) by applying our energy 
principle to this two-fluid model. We show that a rather simple displacement field 
is enough to make the potential energy decrease {6W < 0) and to obtain a tearing 
instability whose growth rate agrees with the asymptotic matching result [6], [7]. Given 
these observations, we extend the energy principle to a nonlinear regime in Sec. HI where 
the displacement (or the magnetic island width) is larger than dg. Without relying on 
perturbation expansion, we directly substitute a form of the displacement map into 
the Lagrangian and attempt to minimize the potential energy W. We show that a 
continuous deformation of magnetic field-lines into a F-shape |25] asymptotically leads 
to a steeper decrease of W than that of the linear regime, which is indeed found to be 
responsible for the acceleration phase. In Sec. O the effect of small dissipation on this 
fast reconnection are considered and implications of our results for sawtooth collapse 
are finally discussed. 



2. Model equations and their Lagrangian description 

We analyze the following vorticity equation and (collisionless) Ohm's law for the velocity 
field V = X V0(x, y, t) and magnetic field B = Vip{x, y, t) x + BqE^: 

^ + [0,VV] + [V^^,^]=O, (1) 

+ [0, - VV] = 0, (2) 

where [/, g] = (V/ x Vg) ■ e^. As noted in Sec. [1], the parameter de denotes the electron 
skin depth, which is much smaller than the system size {d^ <^ L). Since the frozen-in 
flux for Eq. ([2]) is not the magnetic flux ip, but the electron canonical momentum defined 
hy i/je = 'ip — d^V'^ip, the effect of electron inertia permits magnetic reconnection within 
a thin layer (~ d^) despite a lack of resistivity in this model. 

In the same manner as in Ref. [9], we consider a static equilibrium state, 

0(°) = O and 4j'^^\x) = ^pocosax, (3) 

on a doubly-periodic domain D = [—Lx/2, x [—Ly/2, Ly/2] (where a = 27i/Lx), 

and analyze the nonlinear evolution of the tearing mode with wavenumber in the y- 
direction k = 271 / Ly at its early linear stage. For sufficiently small k such that 

TTk^Aa^ = Ll/8Ll < 4 < ^x, (4) 

this instability is similar to the m = 1 kink-tearing mode in tokamaks (which belongs to 



the large-A' regime; see Appendix A). Figure [T] shows contours of ip calculated by direct 
numerical simulation, where e denotes the maximum displacement in the x-direction. 
Our numerical code employs a spectral method in the y-direction with up to 200 modes 
and a finite difference scheme in the x-direction with uniform grid points ~ 10, 000. The 
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Figure 2. Growth of e = e/de with respect to time t = t/ro {de/Lx ~ 0.01 and 
Ly/Lx = 47r) 



growth of e accelerates when e = e/rfg > 1, as shown in figure |2] (which is faster than 
exponential). In accordance with Ref. [H], a strong spike of electric current J = —V'^tp 
develops inside the reconnecting layer and the width of this current spike continues to 
shrink as time progresses, unless a dissipative term is added to ([2]). Therefore, direct 
numerical simulation of ([1]) and (|2]) inevitably terminates when this coherent energy 
cascade reaches the limit of resolution. 

In order to clarify the free energy source of this explosive instability, we solve 
the conservation law ([2]) for ipg = ip — d'^V'^ip by introducing an incompressible fiow 
map Gt : D ^ D, which depends on time and corresponds to the identity map when 
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t = —CO (G_oo = Id). Let {x,y){t) = Gt{xo,yo) be orbits of fluid elements labeled by 
their position (xq, yo) at t = — oo. Then, the velocity field (or </>) is related to Gt as 
dG 

-gf(^o^ yo) = X V0(x, y, t). (5) 

Regarding Gt as an unstable fluid motion emanating from the equilibrium state ([2]), we 
can solve Ohm's law ([2]) by 

^e{x,y,t) = ipe{Gt{xo,yo),t) = ij^^\xo), (6) 

where ipi^^ (x) = (1 + dla'^)'ijjQ cos(Q;a;). Both and ipg (or ip) are thus expressed in terms 
of Gt- By adapting Newcomb's Lagrangian theory [22], we define the Lagrangian for 
the fluid motion Gt as 

L[Gt] = K[Gt] - W[Gt], (7) 

where 

K[Gt] = 11 \Vct>\'d'x, (8) 



2 



' D 



W[Gt] = \ I (I V^|2 + dl\V'ij\') d'x. (9) 
^ Jd 

Then, the variational principle 6 J L[Gt]dt = with respect to 6Gt yields the vorticity 
equation ([1]). 

Since the Hamiltonian corresponds to H = K + W {= const.), we note that W 
plays the role of potential energy and the equilibrium state ([3]) initially stores it as 
free energy. In the same spirit as the energy principle [161 113 , if the potential energy 
decreases {6W < 0) for some displacement map Gt, then such a perturbation will grow 
with the release of free energy. In comparison to the ideal MHD case [22] , the electron's 
kinetic energy (1/2) f^dlJ^d^x appears as a part of the potential energy, because we 
have treated the conservation law of electron's momentum ip^ as a kinematic constraint. 
To avoid confusion, we will refer to this (1/2) d^J'^d'^x as current energy in this work. 

3. Energy principle for linear stability analysis 

In our linear stability analysis, the equilibrium state is perturbed by an infinitesimal 
displacement, Gt{xo, yo) = {xq, yo) + ^(xo, yo, t), where ^ is a divergence-free vector field 
on D. For a given wavenumber k = 2tt / Ly, we seek a linearly unstable tearing mode in 
the form 

sin ky 



^{x,y,t) = V 



emx) 



(10) 



k 

with a growth rate e{t) oc e'^*. We normalize the eigenfunction C,{x) by max \C,{x)\ = 1 so 
that e{t) is equal to the maximum displacement in the x-direction and, hence, measures 
the half width of the magnetic island. 

Upon omitting from equilibrium quantities, V'e°\ J'"'^\ etc., to simplify 

the notation, the eigenvalue problem can be written in the form 

iiVk' + m i'] ' + {^ye + ) i = d'MJ"i + d'Mv'^^^v'm), (11) 
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Figure 3. Test function that mimics the unstable tearing mode 



where should be interpreted as = 9^ — k"^ and the prime (') denotes the x 
derivative. Note, (fTTj) ranks as a fourth order ordinary differential equation (unless 
de = 0) because of the integral operator (1 — rfgV^)"^ on the right hand side. By 
multiplying the both sides of ( |TT|) by ^ and integrating over the domain, we get 
-2/(2) = 1^(2)^ ^Yieie 



-r 



j(2) = 



dx 



v2 



(12) 
(13) 



1 - rf2V2 

The functionals 7^/'^^) and ly(^) are, respectively, related to the kinetic and potential 
energies for the linear perturbation. Hence, by invoking the energy principle [I6] (or 
the Rayleigh-Ritz method), we can search for the most unstable eigenvalue (7 > 0) by 
mmimizmg 

ty(2)//(2) ^ith respect to 1. 
Because we assume the ordering (jlj) that corresponds to the kink-tearing mode, 
the eigenfunction ^ is approximately constant except for thin boundary layers at 
X = 0,±Lx/2 and has discontinuities around them because of the singular property 
of (fTT]) in the limit of {'~f/k),k,de 0. The electron inertia effect would smooth out 
these discontinuities. For this reason, we choose a piecewise-linear test function shown 
in figure [3l In a region containing the boundary layer at x = 0, it is given explicitly by 

{1 for X < -de 

—x/de for -de < X < de (14) 
— 1 for de < X. 

The outer layers at x = ±Lx/2 are equivalent owing to the periodicity and symmetry 
of the problem. By substituting this function into ( |T2l) and ( IT3l) . we can make W^'^'^ 
negative and keep J^^-* finite; i.e., we obtain 



/(2) 



d±^ 



and W^^'> 



9e' 



dgTfj , 



(15) 
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Figure 4. Dependence of linear growtli rate 7 on fc (for d^/Lx = 0.01) 



where Tj^^ 



a^Vo = B' \x=o and we have extracted only the leading-order term (see 



Appendix A for detail). The linear growth rate is therefore estimated as follows: 



7= y/_l^(2)/j{2) = ^ 



776rn 



0.881rn 



(16) 



where Tq"^ = dekrj^^. This result agrees with the general dispersion relation derived by 
asymptotic matching P,[7]. Of course, our analytical estimate of the growth rate depends 
on how good the chosen test function mimics the genuine eigenfunction. Nevertheless, 
the result predicted by the simple function (HM shows satisfactory agreement with the 
numerically calculated growth rate (see figure H]) in the small k region corresponding to 
the ordering (jlj). In the following simulations, we always put kL^ = 0.5. 



4. Variational estimate of explosive nonlinear growth 

Next, we consider the nonlinear phase of the linear instability discussed above. We 
remark in advance that a higher-order perturbation analysis of the Lagrangian (i.e., 
weakly nonlinear analysis) will not be successful, as was already pointed out by 
Rosenbluth et al. for the case of the ideal internal kink mode [H]. For example, the 
Lie-series expansion [26] of ([6]) leads to 

= + ^ ■ + ■ V(C ■ V^(°)) + 0{eydl). (17) 

Thus, such a perturbation expansion easily fails to converge when the displacement e 
(or the island width) reaches the boundary layer width ~ de, due to a steep gradient 
dxi ~ i/de of the eigenfunction inside the layers (see figure [3]). Naive perturbation 
analysis is, therefore, only valid for < e ^ rfg, while e actually exceeds de without 
saturation as in figure [2l 

To avoid difficulties of a rigorous fully-nonlinear analysis, we again take advantage 
of a variational approach. Namely, we devise a trial fluid motion (parameterized by the 
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Figure 5. Deformation of contours of -06 by the displacement map p8 



amplitude e) that tends to decrease the potential energy W as much as possible. When 
such a motion is substituted into the Lagrangian ([7]), it is expected to be nonlinearly 
unstable. 

Owing to the symmetry of the mode pattern, it is enough to discuss the boundary 
layer at x = and, moreover, focus on only the 1st quadrant, < x and < y < Ly/2. 
In a heuristic manner, based on the above linear analysis and simulation results, we 
consider a displacement map Ge : {xo,yo) {x,y), where the displacement in the x 
direction is prescribed by 

' geixo) for (i) < 2/0 < ¥ - ' 



X 



2 



Xo + Tiyo-^] (xo - ge{xo)) for (ii) ^ - I < < ^ ^ - 



2 



/ 

t 2x0 - ge{xo) for (iii) + 1 < < ^. 

The regions (i)-(iii) are indicated in figure [Sl^left) and we furthermore define as 

{e~''xo for < Xo < cie 

4e^"^ for de<xo<de + € (19) 
Xq — e for (ig + e < Xq. 

As illustrated in figure |5l this displacement map deforms the contours of V'e into a 
pattern with F-shaped ends [25]. In a nonlinear regime with c/g <^ e ^ L^,., we find that 
such a deformation decreases the potential energy ([9]) in a manner that is close to the 



steepest descent. Leaving the detailed estimate of 6W to [Appendix B[ our reasoning 
process can be detailed as follows. 

First, in the region (i), the flux ipe of the red area of figure [5](left) is squeezed into 
a thin boundary layer whose width is 2de in figure [5]^right). On the other hand, the flux 
is expanded in the region (iii) and the blue area of figure [5]^left) is almost doubled in 
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Figure 6. Changes of ipe and tp from the equihbrium state ipi — around the 
domains (i) and (iii), due to the displacement map (|18|) with e = Sde- 



figure [5l^ right). The resultant forms of tpe and ip are shown in figure|6l Since the magnetic 
flux ip approximately conforms to tpe except in the neighborhood of the boundary layers, 
both deformations tend to decrease the magnetic energy (1/2) J \ Vip\'^d'^x as when 
de e <^ Lx- This overall loss of magnetic energy is the earmark of collisionless 
magnetic reconnection. 

Inside the boundary layers [i.e., the red regions in figure [SJl^right)], care must be 
taken in representing the formation of the strong current spikes |9], which are observed 
as J = {ipe — 4')/dl in figure EJ^i). These spikes tend to increase the current energy 
(1/2) J d^J'^d'^x in (Q. However, the asymptotic form of the current is approximated 
by a logarithmic function, J ~ r^j^elog \x/de\ for e = e/de ^ 1, and the current energy 
change is, at most, of the second order O(e^). Therefore, in the regions (i) and (iii), 
the dominant contribution of the potential energy decreases at the rate of order O(e^) 
despite the minor increase of the current energy. 

Only in the intermediate region (ii) located between (i) and (iii), does the potential 
energy tend to increase due to the bending of magnetic field-lines over the distance /. 
But, we can minimize this contribution from the region (ii) by taking its width / to 
be sufficiently small: I Ly. We are allowed to use this approximation as far as the 
ordering (jlj) is concerned; Ly is the longest scale length in this ordering and, in fact, 
Lj^ — 7- oo is similar to the behavior of the m = 1 kink-tearing mode. 

By noting that there are, respectively, eight regions that are equivalent to (i) and 
(iii) in the whole domain D, analytical estimates given in Appendix B can be gathered 
into the following: 



SW[G,] ~ SSWii) + 8(5Vr(iii) = -LyTH^dl 



(20) 



for de ^ e Lx- 

To evaluate the nonlinear growth rate of e, it is necessary to estimate the kinetic 
energy. By introducing time-dependence in e{t) via the displacement map f|T8|) . a 
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straightforward analysis (given in Appendix C ) eventually results in 




(21) 



where t = t/rQ. This estimate is not remarkably different from that of the linear regime. 
With these estimates, the Lagrangian ([7]) reduces to 



The equation of motion is, of course, d'^e/dt^ = F(e) with F(e) = —{l/2)dU / de. 

In the linear regime (e <^ de)-, we have already shown that the potential energy 
decreases as U{e) = — 0.776e^ and e(t) grows exponentially. The steeper descent where 
f/(e) = — 0.439e'^ in the nonlinear regime {de <^ e <^ L^) indicates an explosive growth 
of e, namely, it reaches the order of the system size during a finite time ~ Tq. 

We remark that the nonlinear force F{e) ~ O(e^) obtained here is different from 
F{e) ~ O(e^) in the earlier theory by Ottaviani and Porcelli [9]. While similar fluid 
motion around the X and O points is considered in Ref. [9], they directly integrate the 
vorticity equation ([1]) over the quadrant [0, x [0, Ly/2] and arrive at an equation of 
motion d'^e/dv = F{e) ~ O(e^). However, unless the assumed trial motion happens to 
be an exact solution, their treatment may lead to a wrong equation of motion that does 
not satisfy energy conservation. Moreover, their ansatz of the "fixed flow-pattern" is 
also found to be inappropriate in our trial-and-error process. If we try fixing the stream 
function (p throughout the linear and nonlinear regimes as 



as shown in figure [71 With this choice, the potential W does not continue to decrease - 
such a fixed flow-pattern merely circulates the flux tpe from the X point side to the O 
point side via the boundary layer. 

In direct numerical simulation, we have calculated the potential energy U{e) [or, 
equivalently, the kinetic energy {de/diy] as a function of e. As shown in figure [H the 
decrease of U{e) agrees with our scaling and does not support the scaling U ~ — of 



5. Small dissipation 

In this section, we consider the effect of small dissipation by introducing resistivity (rj) 
and electron perpendicular viscosity (/ig) into Ohm's law ([2]); i.e.. 




(22) 



where the normalized potential energy is given by 

U{e) = -{3/TTHog2)e^ + 0{e^) = -0A39e^ + 0{e^). 



(23) 




Ref. [9]. 



dljje 



+ V ■ V^e = -VJ + l^edlV'^J. 



(25) 



dt 
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Figure 8. Potential energy U{e) (where del — 0.01 and LyjL^ = 47r in simulation) 



Both terms on the right hand side only dissipate the potential energy W . For sufficiently 
small ?7 and /ig, we can still employ the energy principle in the manner used to describe 
the resistive wall mode in Ref. [27]. Thus, the energy principle is extended as 

-7^/(2) = ly(2) + lyf + 0(ry^/.D, (26) 



where 



and 



lyf = -- / ' (r/l + iXedW^jf) dx < 0, (27) 



e 
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By substituting the same test function ^ of (fT4|) into W^J, the hnear growth rate ( fT6ll 
is modified to 



where Tg = dl/rj is the electron colhsion time and = dl/i^e is the electron diffusion 
time over a distance de- Small rj and fie, therefore, enhance the linear growth rate. 

Now, let us interpret our result for tokamak parameters. Since Te/ra ~ {pe/deY, 
where pe is the electron gyroradius, the effect of electron viscosity is typically much 
smaller than that of resistivity, Tf,/Td <^ 1, in strongly magnetized plasmas in tokamaks. 

However, the time scales tq and Tg can sometimes be similar in tokamak plasmas. 
For the m = 1 kink-tearing mode in tokamaks, Tq"^ = d^kr^^ corresponds to Tq^ = 
deq'iUAo, where q[ is the derivative of the safety factor q at the q = 1 surface and 
ujao is the toroidal Alfven frequency at the magnetic axis. For sample parameters, 
^Ao = 6.4 X 10^s~^, Tg = 6keV, n = 3.5 x lO^^m"^ and q'l = 2.0m~^, corresponding 
to TFTR experiments that have sawtooth crashes [U [5], we obtain tq = 90/is and 
Tg = 270/is, although the ratio tq/ts = 0.33 can drastically change in proportion to 



By recalling that the resistive layer width 5^ for the case of A' = oo is given 
by Srj ~ {r]/q[ujAoY^^, namely, 5n/de ~ {ro/reY^^, we expect the reconnection to be 
relatively collisionless when tq is shorter than Tg. Indeed, the nonlinear acceleration 
phase is observed numerically for ro/rg < 1. Figure [9] shows instantaneous growth 
rates of e(t) for different values of resistivity ro/rg(oc 77). The linear growth rate 7, 
which emerges at the small amplitude e/rfg = 0.1(-C 1), obeys the dispersion relation 
7^0 = (To/Te + 7To)^^^ obtained by asymptotic matching [2H]. As the amplitude e enters 
into the nonlinear phase e/d^ > 1, acceleration occurs for ro/rg < 1. Since the electron 
skin depth d^ is wider than the resistive layer width 5^ for ro/rg < 1, collisionless 
reconnection governs macroscopic fluid motion. On the contrary, for ro/rg > 1, the 
resistive layer initiates the reconnection process and hence deceleration occurs in figure [H 
which is more like the quasi-equilibrium evolution caused by the resistive kink mode [2]. 

6. Summary 

In this work, we have analytically elucidated the acceleration mechanism for collisionless 
reconnection driven by electron inertia. A variational method based on the Lagrangian 
description of collisionless plasma is shown to be useful especially for predicting nonlinear 
evolution; conventional asymptotic matching does not apply to this problem unless an 
exact nonlinear and unsteady solution is available around the boundary layers. 

We have demonstrated the existence of a nonlinear displacement map that decreases 
the potential energy of the Lagrangian system into the nonlinear regime. No matter how 
small the electron skin depth de, electron inertia enables ideal fluid motion to release 
free energy (~ magnetic energy) of the equilibrium state because the frozen-in flux is 
switched from ip to ipe = 'ip ~ d'lV'^ip, producing a reconnecting layer of width de- In 




(29) 
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Figure 9. Instantaneous growth rates {TQ/e)de/dt versus resistivity to/ts, numerically 
evaluated at several levels of amplitude e/de {de/L^ = 0.01 and Ly/L^ = 47r). 

the large- A' limit, the formation of F-shaped structures connected by a current layer in 
the magnetic configuration is favorable for the steepest descent of the potential energy. 
This descent scales as O(e^) for e = e/de ^ 1 with respect to the displacement e (or the 
island width). The associated explosive growth of e would continue until e reaches the 
system size and leads to an equilibrium collapse during a finite time ~ tq. 

Although our analytical model is too simple to explain all sawtooth physics in 
tokamaks, the time scale of explosion (tq ~ 90/is) that is predicted in this work 
is comparable to experimentally observed sawtooth collapse times [U |5]. However, 
resistivity is not negligible in tokamaks and tends to decelerate the reconnection. Our 
simulations exhibit nonlinear acceleration only for the case of tq/te = r]/dlq[u!Ao < 1, 
which can be fulfilled by the experiments. In more realistic plasmas, the strong current 
spike generated by electron inertia would cause rapid heating of the plasma, which would 
reduce the local resistivity r], and, what is more, would produce runaway electrons. Since 
these effects also act as positive feedback, we expect that sawtooth collapse occurs once 
the acceleration condition ro/re < 1 is satisfied at the q = 1 surface. 

We infer that the state of lowest potential energy is similar to Kadomtsev's fully 
reconnected state (where q at the magnetic axis is go = 1) [I]- But, if dissipation were 
sufficiently small, it would also corresponds to the state of maximum kinetic energy, 
where a strong convective flow remains. As shown in numerical simulations [20l |2T] . 
such a residual flow causes a secondary reconnection and restores a magnetic field similar 
to the original equilibrium (go < !)• If our analytical result is adapted to cylindrical 
geometry, this partial reconnection model will be corroborated theoretically. 

We expect further application of our variational approach to be fruitful for 
describing strongly nonlinear and nonequilibrium dynamics of sawtooth collapses. As 
is well known, finite-Larmor-radius effects and diamagnetic effects would modify the 
island structure and dynamics significantly on the ion scale which is larger than dg. Our 
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approach is feasible even for multiscale problems that require nested boundary layers, 
as long as dissipation is not the dominant factor. Extensions of the present analysis to 
more general two-fluid equations are in progress and will be reported elsewhere. 



Acknowledgments 

The authors would like to thank A. Isayama, M. Furukawa and Z. Yoshida for fruitful 
discussions. This work was supported by a grant-in-aid for scientific research from the 
Japan Society for the Promotion of Science (No. 22740369). P.J.M. was supported by 
U.S. Dept. of Energy Contract # DE-FG05-80ET-53088. 



Appendix A. Linear stability analysis 

In the ideal MHD limit {dg = 0), it is well known that the equilibrium (j3]) has a 
marginally stable eigenmode (7 = 0) which is expressed hj ip = ipo cos K{a\x\ — 7r/2) 
(where k = ^1 — k'^/a^) in terms of'ip = —ip'^. When fc^ < a^, this eigenmode formally 
makes W^"^^ negative; 

^ ^ x=+0 

= _ 2 = -24jlaKsm{K7r) <0, (A.l) 

x=-0 

which also implies that the tearing index is positive, A' = tp' / %()\^Z^q = 2aKtanK| > 0. 
The corresponding ^ = —ip/ip' is, however, discontinuous at x = 0, ±La,/2 and hence 

= 00. It is therefore reasonable to infer that this marginal mode would be 
destabilized by adding electron inertia de <^ L^. 

Note that the integrand of the potential energy f|T3|) is composed of two quadratic 
terms, which are, respectively, positive and negative definite. Since ^ ip' for small 
de, the main role of the electron inertia is to weaken the magnetic tension (equal to the 
former positive term) through the smoothing operator (1 — rf^V^)"^. 

By assuming the ordering @ (in which A' ^ Sa^/vrA;^ is large) and employing the 
test function ^ in figure [HI let us estimate only the leading-order term in ( fT2|) and ( !T3l) . 
For that purpose, we can always use an approximation ^ 9^. Then, ( 1T2|) easily yields 
the estimate of /^^^ in (fTSjl . 

To calculate the potential energy ( |T3l) . we introduce a neighborhood [—do, do] of the 
boundary layer [—de,de], where do is supposed to be a few times larger than de- The 
potential energy in the outer region [—Lx/2 + do, —do] U [do, Lx/2 — do] is estimated by 



='^0 dn 

^ -4^, (A2) 

=-do 



where r^^ = -i/^o^^, because ijj ~ ■j/'o cos(a|x| — 7r/2) in this region. 

Next, we focus on the inner region [—do, do] by using a local coordinate x = x/d^ 
and approximating the equilibrium profile by ^/^g ~ ^Z^' ~ —{de/TH)x. For given 

, -x^ for Ixl < 1 
-\x\ for 1 < X , 
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V 





^(iii) 








































df, de+ e 


de+ 2e 







do - e do (io + e 

Figure Bl. Closeup of three regions (i)-(iii) in figure [5] 



the corresponding is obtained by solving ^jj^ = ^~ d^ip under the boundary condition, 
ijj ^ ipe \x\ oo. This analysis results in 

3 




-(x' + 2) + -e"i(e" + e 



3e-i - e 



\x\ 



-\x\ 



\x\ < 1 
1 < 1x1, 



(A.4) 



where ip{0) ^ indicates that this perturbation causes magnetic reconnection. The 
potential energy inside the layer [—do, do] is calculated as 

1 rdo/de 

-'\^,^p\' + \m 



W, 



(2) 

[-do,do] d 



dx 



do/de 
dr, f I 



H 



.dn 



3 



(A.5) 



where we have neglected e~'^°^'^'' by making do larger than de to some extent. Since other 
boundary layers at x = ±Lx/2 can be treated equivalently, the total potential energy 
on the whole domain is estimated as (IT^ . 



Appendix B. Estimate of potential energy change 



Here we estimate change of the potential energy that is caused by the nonlinear 

displacement map given in (1181) . In a way similar to that of the linear analysis (see 
Appendix A ), we introduce a domain [0, c/q] x [0,Ly/2] where do is now taken to be 
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somewhat larger than 2de + 2e. This domain is deformed by the map as shown in 
figure IBlt where we refer to shrinking and expanding domains as 



D 



(iii) 



[0,do 
[0,do 



e X 



+ e] X 



'4 2 

hi A_ L hi 

4 ^ 2' 2 



(B.l) 
(B.2) 



respectively, and their intermediate domain as -D(ii). The potential energy will not 
change significantly outside of these domains, because the fluid is simply subject to 
parallel translation along the x direction. Moreover, we are only interested in the 
domains Z)(i) and -D(iii), on which efficient decrease of the potential energy is observed 
as follows. 



Appendix B.l. Potential energy on 

The equilibrium is approximated by 'ilji^\x) ^ tp'^^^x) ~ V'o(l ~ a^x^/2) and it is 
deformed into ipe{x,y,e) = ipi^\g~^ (x)) on Therefore, 6ipe = "^e — 'ipe'^ is given by 

r 1 

2' 



1 



:i - e'nx^ 



for < a; < dpC ^ 



X 



rfe I log — + 1 I + e 

dp 



-ex 



2 2 

+ — for (ieC"^ < X < de 



for dp < X. 



(B.3) 



For large e = e/de ^ 1, we can neglect the innermost region < x < dpC ^ and the 
asymptotic form of Sipp contains a logarithmic function as follows. 



5Ve 



1 x^ 
'1 + e) logx - -(1 + e)^ + — forO<x<l 



(B.4) 



for 1 < X, 



where x = x/cig, and the corresponding spike is recognized in figure |6]|^i) . By solving 
(1 — dDdtp = 6ipe for 6ip, the change of current 6 J = —dl{6ip) turns out to be 

1 f -(e + l)Ec(x) - 1 + ci cosh(x) for < x < 1 
6J= — {_.^^' ^ ' (B.5) 

th \c2e for 1 < x, 

where we have defined Ec(x) = [e'^Ei(— x) + e~^Ei(x)]/2 using the exponential integral 
Ei(x) = p.v. J^^i^e"^ /s)ds. The coefficients Ci and C2 are matching data at the interface 
X = 1, which are linear functions of e as follows; 

(e + l)[Ee'(l) + E,(l)] + l 



(B.6) 
(B.7) 



C2(e) = (e + 1) [cosh(l)Ec'(l) - sinh(l)Ec(l)] - sinh(l). 

Since Ec(x) ^ log |x| near x = 0, a strong current spike develops in the form of 
logarithmic function, as noted earlier in Ref. [9]. However, the asymptotic form of 
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6J = J — J^") remains square-integrable and, hence, the current energy change in is, 
at most, of the second order; 



do 



rf2(|J|2_|j(0)|2)^^ 



(B.8) 



H 



On the other hand, the magnetic flux = ijj^'^'^ + 5ip is free from such logarithmic 
singularity and its derivative. 



for < a; < 1 
for 1 < X, 



(B.9) 



-(1 + e)dx [logx — Ec(x)] — Ci sinh(x) 
-(x + e) + C2e"^ 

is again square-integrable and linearly depends on e. Hence, the leading-order estimate 
of magnetic energy is simply 







do-e j3 

Idx^Pl'dx = ^ 

H 



T 



(do/de)—i 



(x + efdx + 0(e^ 



'H 



• j3 3 

f - ^ + Oie^dl) 



(B.IO) 



which decreases as for e ^ 1. This decrease of the magnetic energy dominates the 
increase of current energy (IB.Sp . Therefore, the potential energy change on is found 
to be 



hi 
4 



'H 



6 



(B.ll) 



Appendix B.2. Potential energy on -D(m) 



For the purpose of estimating 5W on -D(iii) to leading order, one may approximate the 
inverse map of fllSp as 

X 

2 (B.12) 



Xq 



X 



for < X < 2e 
for 2e < X, 



for large e. The equilibrium flux il)f\x) ~ '?/'o(l ~ is expanded by this outflow 

and is deformed into a flat-topped shape [see flgure [6t^iii)]. In the same manner as for 
the domain we flrst obtain 



th 1 2ex 



for < X < 2e 
for 2e < X, 



(B.13) 



and calculate the current change as follows: 



6J = — 



3 
'i 



2 4/2^ 



! 3~ 
2^4/2 



x-2e 



for < X < 2e 
for 2e < X . 



(B.14) 
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By keeping the smallness of e ^ in mind, we confirm that the current energy change is 
again of the second order O(e^). The asymptotic form of d^ip is estimated by 

X I e 3\ 1 



4^ 



■ + - I -e' 
d. \ 4 V 2 4 / 2 



x-2e 



th 



[X 



3\ 1 



-x+2e 



2 4/2 

and the magnetic energy is also found to decrease as e^; 



for < a; < 2e 
for 2e < X, 



do+e 



\dx'ijj\'^dx 



4_ 

2rl 



lb 



{x-efdx^O{e^ 



2e 



6 12 ^ 



(B.15) 



(B.16) 



The flat-topped region of '?/'e(— V^) corresponds to the magnetic island, on which the 
magnitude of d^'^l^ obviously decreases. The potential energy on -D(iii) therefore decreases 
as follows: 



(iii) 









■ e3 




2y 


4 


_~12 



?2^3n 



(B.17) 



Appendix C. Estimate of kinetic energy 

Here we estimate the kinetic energy K\G^(t)\ of the displacement map Ge(f) given in 
(IT5]) . where only e{t) is assumed to be time-dependent. By invoking figure [BT] again, the 
dominant part of kinetic energy turns out to exist in the domains -D(i), -D(ii) and -D(iii). 
Since is ignored in this work (by assuming / ^ Ly), we exhibit only the results for 
il'(i) and -D(iii) as follows. 



Appendix C.l. Kinetic energy on 

Owing to our special choice of g^^, the x-component of the velocity field on /^(i) is simply 
given by 



X 



de 



— ioT < X < dp 



(C.l) 



-1 for de < X. 

By solving the incompressibility condition d^v^ + dyVy = under appropriate boundary 
conditions, the ^/-component of the velocity field is found to be 

y 



Vy{x,y,e) 




for < X < de 
for dp < X. 



(C.2) 



This Vy dominantly contributes to the kinetic energy on -D(i), which is readily estimated 
by 
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Appendix C.2. Kinetic energy on -D(iii) 

We can go through the same procedures as for but the analysis is somewhat 

comphcated by the fact that the inverse map x (-> should be dealt with as an implicit 
function. In terms of the unperturbed position xq, the velocity field is expressed by 

' for < xo < de 



de dx 
dt de 




Vy{x,y,e) 



de 1 



2 - e-^ 

gXQ-e-l 





for de < Xq < de + e 
for de + e < xq, 

for Q < Xq < de 



(C.4) 



for de < Xq < de 
for de + e < Xq. 



XC.5) 



Using the change of variables from x to Xq, the kinetic energy on -D(iii) is therefore 
estimated as 



(iii) 



dy 



do ]^ 

dxo-{vl 



21og2 - 1 fL 



Qde 



4 



n, dx 
) 

^ dxo 

2 



(C.6) 



where we have neglected e for large e ^ 1. 
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